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ABSTRACT 


The M-Layer program tracks the constant phase lines Im{D(q)} = 0 and looks for 
their intersections with the lines Re{D(q)} = 0 for the locations of the zeros of the mode 
function D(q). These two types of constant phase lines are tracked and plotted over a 
search region which contains modes having a range attenuation rate of no more than 5 dB 
per km. Several new parameters for use in mode search are deduced from the results and 
some old ones are verified. Future studies in waveguide mode propagation theory 
pertaining to atmospheric ducts may benefit from this work. An improved mode search 


strategy is also proposed. 


TABLE OF CONTENTS 

I. INTRODUCTION (0°) BOUT. oe tee ae eter cee ce, ne 1 
A. THE WAVEGUIDE MODE THEORY OF PROPAGATION ....... 2 

B. MODE SEARCH PRINGIPRE . 2255... FR eee ‘i 

II. MODE FUNCTION IN THE COMPLEX q,j/t; PLANE................ ul 
A, “MODE*IEOCATIONS oo... eo Oe ce 0 + 11 

B. PHASE LINE TRACKING. Sie oi-...-..++...00 13 

1. Downward Traekiig ..............,....... Ie 

2. Upward Tracking fia... dame -..........:.. rr 18 

Il, ANALYSIS AND CONCLUSIONS¥e= wo 5005s...) 21 
A. MODE SEARCH PARAMETERS ................00 cee eees 21 

B. MODE SEARCHSTRAGEGY 24a. . (2. ....... eee 24 

1, Track Pemminatioa,. . 9 cs ee ss 6 «oes. oi 24 

2. Search Temmniionge: -- eer ueree......2 cee 25 

3. Track Duplication Avoidance ...........+.......: eee 27 
APPENDIX A: MODE LOCATIONS IN THE COMPLEX q,j/t; PLANE ..... 28 


APPENDIX B: SUBROUTINE FNDMOD AND FZEROX ................ 49 


APPENDIX C: CONSTANT PHASE LINE TRACKING FLOW CHARTS ..... 61 


APPENDIX D: CONSTANT PHASE LINES IN THE COMPLEX q, ,/ty 


fee NDIX E; SEPARATION OF Im{D(q)}=0 LINES ................... 75 


APPENDIX F: DIFFERENCE BETWEEN CONSECUTIVE Re(q,/t;) VALUES 84 


ENC SV RIS, SUE oe a 93 


ACKNOWLEDGEMENTS 
I would like to express my sincere gratitude to my advisor Dr. Hung-Mou Lee for 
his professional guidance, assistance and support throughout the thesis preparation period. 
I thank my second reader Dr. Lawrence J. Ziomek for his careful review and constructive 
suggestions to the thesis work. 
Special thanks to my wife Mei- Ying Huang and my son Yu-Heng Che for their love 


and unending support which allowed me to successfully complete my graduate education. 


Vi 


I. INTRODUCTION 


The M-Layer program developed by NOSC (Naval Ocean System Center) was 
documented by Yeoh [Ref. 1] at NPS. It was later extensively revised by Lee and Han 
[Ref. 2] for improved accuracy and speed. The new FORTRAN code is now identified 
as the NPS version [Ref. 3] under the auspices of NRaD (Naval Command, Control and 
Ocean Surveillance Center, RDT&E Division), the current name for NOSC. In this 
version, the logical structures and numerical algorithms for following the constant phase 
lines of the mode function and for computing the mode locations were almost completely 
‘rewritten. The overall plan to first partition the search region into rectangular areas and 
then circulate around these rectangles to search for the modes, i.e., the mode search 
protocol, was left unchanged. Lee and Han [Ref. 2] recommended that a more direct 
mode search protocol should be devised to locate the modes more expediently. If possible, 
such an approach should locate the modes in the order of ascending range attenuation 
rates. 

To design a more efficient mode search protocol, a better understanding of the mode 
function beyond the known locations of the modes is required. In this thesis, the 
analytical property of the mode function is investigated. Programs are written to track and 
plot the constant phase lines along which the mode function is either real or purely 


imaginary. From the results, a new strategy to locate the modes is recommended. 


In the remainder of this chapter, the background theory and some of the notations 
used in this thesis are introduced. The basic theory presented in Sections A and B follows 
Ref. 3 closely. The results of this research are presented in Chapter II. In Chapter III, a 
new mode search protocol is proposed based on the findings of this thesis. The relevant 


parameters for use in this new approach are also discussed. 


A. THE WAVEGUIDE MODE THEORY OF PROPAGATION 

Trapping of electromagnetic (EM) waves in the modes supported by a duct is the 
dominating factor in over-the-horizon propagation. The computer program M-Layer 
searches for these modes and computes the electric field from the Hertz vector. Assume 


a vertical electric dipole 


J = 4n28(x)8(y)8(z-z,) 


at a height Z=Zy above the surface of the earth. Following Freehafer [Ref. 4], the Hertz 
vector of the EM fields can be written, in the cylindrical coordinates (p, $, z), aS a sum 


of contributions from individual modes [Ref. 1]: 
II(p,z)=-nj), Hy (B,P)8q(Z)Bq(2) (1) 


where 8, is the wavenumber of the m-th mode and is independent of the coordinates; 
Zr is the height of the transmitter above the ground, Hy is the Hankel function of the 


second kind, which represents an outgoing wave in the radial direction when the time 


dependence Jt js adopted, and g_, is the height-gain function of the m-th mode, 
normalized so that the integral of its square over all heights equals unity when either an 
electric or a magnetic dipole source is used. The height-gain function satisfies the 


equation 
ae { k?m (z) p Zé, =0 ( ) 


where k is the free-space wavenumber, m(z) is the modified index of refraction, and 


m*(z) 1S approximated with a continuous, piecewise linear profile having I layers: 


m*(z)=m; +0,(z-z,) qa, “itstell. (3) 


In practice, m(z) deviates only slightly from unity in the troposphere. The modified 
refractivity M(z) = [m(z)-1]x109, as shown in Fig. 1, is commonly used. The slope of 
M(z) is approximately ax10°/2. 


In the i-th layer, the height-gain function is given by 


g,,(2)=B(B, IC (BK (4p) + (Gq) Z)S$2SZ,,, (4) 


where q,,; 1S a dimensionless linear function of height z with k, Baw m;, and @; as 


parameters: 


ee 
layeri f 


minimum M(z) 
x bala od 


piecewise 


linear 
aitict segments 
height 
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Figure 1. Typical modified refractivity M(z) profile. 





|(&Y | 2 B, (5) 
Qmt = fs Le a ary 
The functions k,(q,,;) and k>(q,,;) are given in terms of the Airy function Ai: 
ky Qn) =J2(12)"*Ai(-q,,€7") (6) 
and 
k,(q,,)=2(12)""Ai(-q,,). (7) 


For z> 2,1, m?(z) is extended continuously upward at a constant slope 
commensurate with the effective radius of the earth. This slope equals 2.36x107/ (m7!) 
for the four-thirds effective earth radius model. Thus, the height-gain function g,_, is again 
given by Eqs. (4) through (7) for z > z;,), with the constant Cy,; set to el4t/3 It 
represent an upward going wave as z becomes large. Beneath the "flattened" earth surface 
where z < z, =0, m?(z) is deduced from the dielectric constant and the conductivity of 
sea water, which are set to 80.8869 and 4.64 (Si/m) respectively. Here the Hertz vector 
is represented by a downward propagating plane wave in this region. 

The conditions that the tangential components of the electric and magnetic fields 
are continuous across the layer boundaries determine the coefficients C;. Starting from 
the top layer down, or "integration-down" [Ref. 1], every C; is uniquely determined by 
C41 at z = Z,1- Thus, a set of all C; coefficients is determined without considering the 


boundary conditions at z = z,. On the other hand, starting from the lowest boundary at 


Z = Z, up to Z = Zy, or “integration-up", a second set of C; coefficients is determined 
without applying the boundary conditions at z = zy,;. That these two sets of C; 
coefficients are identical is a manifestation that Ban is the wavenumber of a mode. 
Mathematically, this consistency condition, sometimes called the guidance condition 
[Ref. 5], can be expressed as the vanishing of a determinant D called the mode function 
[Ref. 6]. The elements of this determinant consist of linear combinations of ki Gine and 
k5(q,y)) and their derivatives with respect to height at the boundaries. The M-Layer 
program searches for the zeros of the mode function to find the wavenumbers. Once a 
wavenumber is obtained, the height-gain function of the corresponding mode can be 
computed. Note that the normalization condition on g,, and the boundary conditions, 
together with the C; coefficients, determine the B; coefficients to within an overall sign. 
This sign need not be resolved because only products in the form of g,,(z)g,,(z7) from 
each mode contribute to the Hertz vector. 

Since the mode function D depends on the wavenumber B_, only through the values 
of q,,j at the layer boundaries, it is more convenient to consider D as a function of the 


variable q given by 


3 2 bs 
q = a3 mee F (8) 

a, ke 
and to search for the zeros of D(q) in the complex q plane. The m-th zero is designated 


as q,,, and is called a q-eigenvalue and B an is then deduced from q,,, by inverting Eq. (8) 


for B. In the NPS version, q,,, is ordered in ascending attenuation rate of the mode, which 


is approximately proportional to the imaginary part of B,_. 


B. MODE SEARCH PRINCIPLE 

M-Layer searches for modes which have attenuation rates below a specified value. 
At ground ranges more than several wavelengths away, this attenuation rate is 
proportional to the imaginary part of the wavenumber B_,, as can be seen from Eq. (1). 


In the upper complex q plane, for values of q such that 


<<1 (9) 














m4—B/k is approximately proportional to q because the absolute value of m, and, hence, 
that of B/k, are both nearly unity. Thus, the limit on attenuation rate imposed on the 
imaginary part of B can be translated approximately into an upper bound! on the 
imaginary part of q which defines a strip in the complex q plane called the search region. 
This region is covered with layers of identical square meshes whose sides are parallel to 


the imaginary q axis and have a length equivalent to an attenuation rate of 1/32 dB per 


1 Since the absorption is small in air, the modified index of refraction m(z), and 
hence @., are considered as real quantities in the following discussions. In the actual 
FORTRAN code, they are declared as complex variables. 


2 This is the default value of the NPS version which can be adjusted through editing 
the variable "dmesh" in an ASCII input file. 


kilometer. The lower edges of the meshes in the lowest layer sit along the real g axis? 


The meshes in the layer just below the top one contain the upper bound of the search 
region, with the top layer providing some allowance for numerical inaccuracy. The 
program tests each mesh square for the presence of zeros of the mode function D(q). 
To search through all the meshes, the program first divides the mesh covering the 
search region into "contour rectangles" with equally spaced vertical lines parallel to the 
imaginary q axis. These vertical lines contain the edges of stacked square meshes and are 
separated by a distance 160 times the side of a mesh. The search commences at the top 
left corner and moves counterclockwise around each "contour rectangle" and begins with 


the one whose left edge is defined by 


mot)” 


where m,,,;,, 18 the minimum value of the modified index of refraction profile and d,.,i 


= Re (m} -majq) ~2x10°6 Re(M, -M,in) =2*10°°d min a) 





n 


is shown in Fig. 1. The justification for this choice to begin the mode search is as 
follows: From Eq. (5) and optics, for a wave to propagate freely in space, the dispersion 


relation of the Maxwell equations requires that k*m(z) > p assuming that both 


m’ 


quantities are real. Thus, the smallest so to support a trapped wave will occur when 


3 The NOSC version extends the lower edge slightly below the real q-axis. This 
causes problems in some situations [Ref. 2]. 


pe just exceeds the minimum of m(z)*. It is not surprising that an argument based 
on optics works within a waveguide mode formulation which is low frequency in 
principle. Lee [Ref. 7] has demonstrated that the earth-flattening approximation actually 
links Mie’s low frequency oriented spherical harmonics to Fock’s high frequency 
diffraction theory through the uniform asymptotics. The mere introduction of an unlimited 
ground range into the Maxwell differential equations incorporates the global feature of 
the radius of the earth into the local equations. 

After the search over the initial rectangle is completed, the program goes on to 
search the neighboring rectangle to the left. If a specified number” of consecutive 
rectangles of decreasing real q values have been searched without turning up any zero of 
D(q), the program changes direction and starts to search the rectangles to the right of the 
initial “contour rectangle" one by one, with increasing real part of q. After failing 
consecutively to locate any q-eigenvalue again over the specified number of rectangles, 
the program assumes that no more zeros of D(q) exists in the search region. The mode 
search is considered complete and the procedure is terminated. If the array for storing the 
q-eigenvalues is filled up before the search is completed, the search is terminated with 


an error message. 


4 Note that assuming the same real part, an imaginary component of B,, reduces the 
real part of B“ , which may enable the wave to propagate in the z-direction. 


> This number of consecutive rectangles is an adjustable input variable "nstop" in the 
NPS version. The default value is 2. 


The search for zeros of D(q) makes use of the fact that a real valued function 
changes sign when it crosses a simple zero. Since a zero of the complex valued function 
D(q) is where both its real part and imaginary part vanish, a necessary condition for a 
point q,,, to be a zero is that it is the intersection of two curves defined by Im{D(q)} = 0 
and Re{D(q)} = 0. M-Layer moves along each side of a "contour rectangle” while 
searching for a sign change in Im{D(q)} across an edge of a mesh bordering the side of 
this rectangle to determine that a line of Im{D(q)} = 0 has been encountered. The search 
then follows this line into the meshes within the "contour rectangle", checking each mesh 
to see if a curve Re{D(q)} = 0 enters the mesh being inspected. All these steps make use 
of the assumption that the zeros of D(q) are simple. Once both the curve Im{D(q)} = 0 
and the curve Re{D(q)} = 0 are found to be present within a mesh, the locations of their 
possible intersections are estimated. 

The rule for sign change becomes inconclusive if a zero of Im{D(q)} = 0 or if 
Re{D(q)} = 0 happens to fall on a corner of a mesh, and a remedy is required. In the 
NPS version, whenever the real part or the imaginary part of D(q) vanishes on a corner 
of a mesh, the phase angle of D(q) is rotated by 2752 radians. This maneuver effectively 
shifts q by a small amount to resolve the sign ambiguity. 

To estimate the locations of zeros of D(q) within a mesh, D(q) is assumed to be 
well approximated in this mesh by its four-term Taylor series expansion. The unshifted 
values of D(q) at the mesh corners determine this cubic polynomial uniquely. Cardan’s 
formulas [Ref. 8] are used to locate the zeros of this polynomial before retaining only 


those lying within the mesh. 
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Il. MODE FUNCTION IN THE COMPLEX q, ,/t,; PLANE 


A. MODE LOCATIONS 

M-Layer searches for the zeros of the mode function D(q) in the complex q plane 
by tracking the constant phase lines of D(q) along which the mode function is either real 
or purely imaginary. The zeros found are called the q-eigenvalues as defined by Eq. (8) 
and are denoted as q,,. These eigenvalues are saved in an ASCII file and are utilized later 
for height-gain function computations. In order to design a strategy for locating these 
zeros, the mode locations are plotted. 

It is clear from Eq. (8) that the variable q varies with 1, the slope of m(z) in the 
lowest (first) layer. Since a, depends on how the continuous, piecewise linear 
approximation to m2(z) is made, while both m, and B.,, are, in principle, dependent only 
on the actual profile, q will vary strongly with @; and fail to provide information on the 
wavenumber of the mode directly. A more suitable variable to use is, from Eq. (8): 


fz) ont-E. (11) 
aA 


Since q is given the variable name of q,, and (k/a Je! 3 is given the variable name of t, 
in the M-Layer FORTRAN code, this variable is identified as 4, /t, throughout this 


thesis. Since m, is real and both it and B/k deviate from unity only slightly for all cases 


11 


investigated, q, ;/t; is approximately 2(m, — A/k). Thus, its imaginary part is directly 
proportional to —-f, the attenuation rate of the mode. Furthermore, q, /t; does not depend 
explicitly on a. Removing the t, factor from q,, makes it possible to compare results 
from different evaporation duct profiles meaningfully. 

Plots of the q-eigenvalues as individual points in the complex q,j/t, plane are 
included in Appendix A. For the 20 meter duct, which will be used as the representative 
case for discussions in this thesis, a line linking all the eigenvalues in the order of 


increasing attenuation rate is shown in Fig. 2. Because of the long excursion between 





Figure 2. q-eigenvalues in the qj, /t; plane, connected in the order of ascending 
attenuation rate (20 m duct). 
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many pairs of consecutive modes in the upper part of this figure, it appears that the 
intuitive approach of starting with the mode of the lowest attenuation and continually 
looking for the mode with the next higher attenuation rate may not be the most efficient 
strategy. Furthermore, there seems to be a fork near the middle of the figure which 
suggests that there can be more than one constant phase line passing through all the zeros 
of the mode function. These problems suggest that, in order to design an efficient mode 
search strategy, the behavior of the mode function in the complex q, ;/t, plane has to be 


investigated more thoroughly. 


B. PHASE LINE TRACKING 

The M-Layer program follows a line along which Im{D(q)} =0 until a zero is 
encountered at its intersection with another line along which Re{D(q)} = 0. Since there 
is no way to determine the phase of the constant phase line linking two adjacent zeros a 
priori, following these two types of somewhat arbitrary but easy to compute phase lines 
which have constant phases of 0 or x, and 2/2 or 31/2, respectively, are still the best way 
to program a computer to locate the zeros systematically. It is thus necessary to find out 
the actual distribution of these constant phase lines in the complex q, /t, plane. 

1. Downward Tracking 

In the original design, along the real q axis, M-Layer divides the search region 

into "contour rectangles," each of which spans 160 meshes horizontally. From the 
observed locations of the modes, it is found desirable to track and plot the constant phase 


lines along Im{D(q)} =0 and Re{D(q)} =0 over a little more than six "contour 


1S 


rectangles", roughly three to the right and three to the left of the starting real q coordinate 
given by Eq. (10), which is a variable called q,.., in M-Layer. The subroutine FNDMOD 
and FZEROX are modified to track these constant phase lines. A listing of these modified 
routines for tracking the lines Im{D(q)} = 0 beginning from the top edge and moving 
downwards into the search region is included in Appendix B. Their flow charts are given 
in Appendix C. Only minor changes are required to track the type of lines Re{D(q)} = 0, 
or to track these lines from the bottom of the search region upwards. The program first 
searches for a sign change in Im{D(q)} or in Re{D(q)} along the top edge of the search 
region, starting at the coordinate -512 mesh sizes from q,,.,, until a distance of 1024 
mesh sizes is covered. When a sign change is observed, a constant phase line is 
recognized as passing through the particular mesh square and the program writes the 
complex q value of the lower-left corner of this mesh square into an ASCII file identified 
by whether the mode function is real or imaginary along the line, and the order this line 
is found. The program then moves from the top edge downwards to follow this constant 
phase line until it exits the search region. The q values of the lower-left corner of the 
mesh squares along this phase line are also written into the same ASCII file to be plotted 
later using MATLAB/386. The constant phase lines for the mode functions of the 2, 4, 
6, 8, 10, 20, 30 and 40 meter evaporation ducts are obtained and plotted. They are 
included in Appendix D. Tracking and plotting these constant phase lines are extremely 
time consuming. The CPU time required to track the desired constant phase lines for each 
duct using an Intel 80486 based PC running at a 33 MHz clock rate is listed in Table 1. 


The time required to plot those lines for each duct using an Intel 80386 based PC running 
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at a 16 MHz clock rate is also listed. 


TABLE 1 
CPU Time for Tracking and Plotting Constant Phase Lines 
“racking racking plotting 
Im{D(q)} = Re{D(q)} = both lines 


(hr:min:sec) (hr:min:sec) (hr:min:sec) 


——— 9:08:14 9:06:45 0:36:49 
ei 8:39:04 8:37:54 0:36:53 
45:07:33 45:04:09 4:01:29 





The constant phase line plot for the mode function of the 20 m evaporation 
duct is given in Fig. 3. The solid lines represent those having Im{D(q)} = 0; the dotted 
lines represent those having Re{D(q)} = 0. Every intersection of these two types of phase 
lines is the location of a q-eigenvalue scaled by tj, thus indicating the presence of a 
mode. Note that every solid line intersects with a dotted line except for many of those 
running from the top through the bottom of the search region. None of the solid lines 
intersect with other solid lines, neither do the dotted lines. This is seen clearly in Fig. 4, 


which magnifies the lower center portion of Fig. 3. One feature that can be recognized 
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the q, ;/t; plane (20 m duct). 


lines in 


Figure 3. Constant phase 
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Constant phase lines in the ql1/tl plane for 20 m duct 





Figure 4. Magnified lower center portion of Fig. 3 showing no intersection of phase lines 
of the same type. 
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immediately in all the plots in Appendix D is that the starting real q coordinate given by 
Eq. (10), which is marked by an "x" along the top edge of each plot, correlates very 


closely with the mode of minimum attenuation. 


2. Upward Tracking 

One mode of low attenuation rate is missing from Fig. 2 compared to those 
found in Ref. 2. Several of them are also missing from those cases of greater duct heights 
plotted in Appendix D. It is evident that tracking the constant phase lines from the real 
qj /t; axis upwards is necessary. For the 10 meter and the 20 meter ducts, Fig. 5 shows 
several constant phase lines starting out of the lower limit of the search region and 
returning to the lower-half complex q, ;/t; plane. Fig. 6 shows the presence of such lines 
for the 30 meter and the 40 meter ducts. These constant phase lines have not been 


observed in those cases of lower duct heights. 


18 


c0 





Figure 5. Upward going constant phase lines in the q,,/t plane for the 10 meter (top) 


and 20 meter (bottom) ducts. 
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11/t) Plane for the 30 meter (top) 


Figure 6. Upward going constant phase lines in the q 


and 40 meter (bottom) ducts. 
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Ii. ANALYSIS AND CONCLUSIONS 


A. MODE SEARCH PARAMETERS 

The implementation of a mode search procedure requires several parameters. The 
search region has to be defined first. In M-Layer, the desired maximum range attenuation 
rate of the modes which support wave propagation is treated as an input parameter called 
Ajogsr Siven conventionally in dB per km. This parameter determines the region over 
which the modes are to be searched as explained in Chapter I. In the complex q plane, 
under the assumption that both m, and B/k of interest are close to unity, the search region 
is bounded from above by the FORTRAN variable Atop given by 


oN 
Trop ~ kxlog ioe 


(12) 


Qioss 4 


The location at which to start the mode search, test» 1S another parameter 


determined by M-Layer. Based on optical considerations, the M-deficiency, d of 


mmin 
Fig. 1, which is the amount of decrease of the modified refractivity from its value on the 
sea surface to its minimum value in the air, is a measure of the capability of the duct to 
trap electro-magnetic waves. As explained also in Chapter I, this quantity determines the 


location for the program to start searching for modes. The validity of the optical argument 


and the usefulness of Eq. (10) have been proved by this work, as pointed out in Chapter 
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To search for zeros of the mode function, M-Layer first covers the search region 
with identical mesh squares. It then follows the lines of the type Im{D(q)} = 0 through 
each mesh square, checking for indications that a curve Re{D(q)} =0 is entering the 
same mesh so that an intersection is possible. The term "mesh size" will denote the size 
of the edges of the mesh squares. It is the size of the step taken by the program to 
advance through the search region. It also determines the initial resolution of the location 
of a mode. The choice of the mesh size should strike a balance between the desire for a 
speedy completion of the search process and the requirement in mode locating accuracy. 
As reported in Ref. 2, for all the evaporation ducts considered, the choice of a mesh size 
equivalent to an attenuation rate of 1/32 dB per km appears to be optimal, that is, all 
modes can be found for all cases investigated in Ref. 2 without experiencing 
extraordinarily long computation time. Consider this mesh size as the default and call it 
qq in the complex q plane. Under the same assumption for deducing Eq. (12), the value 
of q/t, at 9.6 GHz equals 3.58x107°, Figure 7 shows the separation between two 
adjacent constant phase lines of the type Im{D(q)} = 0, indicated in the figure as solid 
lines, along the top edge of the search region for the 20 meter duct. The minimum of 
these separations in terms of qj ;/t, is about 2.81077, which corresponds to a spacing 
of a little less than eight mesh squares apart. The minimum separation between adjacent 
Im{D(q)} = 0 and Re{D(q)} = 0 lines is thus about four meshes. Constant phase line 


separations for all evaporation ducts considered are included in Appendix E. Note that the 


minimum separation is almost constant for ducts higher than 20 meters. It increases 


22 


ay 
Oo 
=) 
so) 
= 
>) 
N 
ww 
- 
oO 
2 
=| 
= 
Zz. 
© 
= 
a 
oO 
7) 
SS 
<= 
Qn 


<5 Ls a a 
ao g (oe) Cr) Cr) (oe) 
(11/1 [b)ay jo souarajsIG 


Figure 7. Separation of Im{D(q)} = 0 lines along the top edge of the search region 
(20 m duct). 
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slightly as duct height is decreased. 

As the program follows an Im{D(q)} = 0 line into the search region, a neighboring 
Re{D(q)} = 0 line moves closer and enters into the same mesh square. This causes the 
program to invoke the root finding routine to determine the possible locations of zeros 
of the mode function within this mesh square. If the mesh size is too large and the 
Re{D(q)} = 0 line enters the mesh square long before the actual intersection takes place, 
the root finding routine will be prematurely invoked many times and the probability of 
producing false modes is increased. This explains why a larger mesh size sometimes will 
increase the execution time. 

Given these three parameters Atop’ Atest’ and qq, a mode search strategy which does 


without the "contour rectangles” is proposed. It should improve the efficiency of M-Layer. 


B. MODE SEARCH STRATEGY 

From the constant phase line plots of Appendix D, a strategy to search for the mode 
can be drawn: move along the top edge starting at q,..;. Search first toward either the left 
or the right, then reverse course to search along the other direction. After the search along 
the top edge is completed, search the lower edge from one end to the other. In what 


follows, the implementation of this strategy will be discussed. 


1. Track Termination 
The tracking of an Im{D(q)} = 0 line naturally ends if the top edge or the 
bottom edge of the search region is reached. On the other hand, the search region is 


unbounded to the right and left of rest: It 18 Convenient to retain the feature in the 
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original program to set a limit on the number of steps allowed to follow a constant phase 
line. From Fig. 3, this limit can be set to 2.5 times the number of mesh squares between 


the vertical limits of the search region. This number equals 2.5x32xa),.... 


2. Search Termination 

M-Layer has to determine that no more modes within the search region is to 
be found and to terminate the search for modes. When searching along the top edge of 
the search region for the constant phase lines Im{D(q)} = 0, the separation between the 
real parts of the mode eigenvalues is of interest. Figure 8 shows the separation in the real 
part of neighboring q-eigenvalues scaled by t), plotted against their locations q,,/t, along 
the real qj /t; axis for the case of the 20 m duct. Similar plots for all duct heights are 
grouped in Appendix F. Excluding the lowest point which involves the mode obtained via 
searching along the lower edge, the separation between neighboring modes is erratic with 
an upward trend away from the center of the figure, which is close to the search starting 
position. 

The increase in distance between two modes towards the end of the range 
searched makes it difficult to implement an adaptive mechanism to terminate the search. 
On the other hand, Fig. 3 shows that almost all phase lines entering the search region 
from the top that contain a mode within this region are bunched together. Therefore, the 


parameter n set equal to four and may be adjusted, can be used to stop the search 


stop” 
along the top edge when four consecutive constant phase lines tracked turn up no mode. 
The search along the real q axis can be confined to within the end points of 


the search along the top edge. 
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Figure 8. Difference between consecutive Re(q, /t)) values (20 m duct). 
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3. Track Duplication Avoidance 

When the program searches step by step along the top or bottom edges for 
sign changes in Im{D(q)} to start tracking the constant phase line, the exit of a constant 
phase line from a previous track will be encountered. Entering into the search region at 
such a location will simply retrace a constant phase line which has already been examined 
for the existence of a mode. Thus, the exit locations of the constant phase lines must be 
recorded and checked whenever a sign change in Im{D(q)} is found before deciding to 
follow the constant phase line into the search region. Along the top edge, a record of four 
updated most recent exit locations from the top edge should be adequate. The locations 
of exits from the bottom should be kept for use during the search along the lower edge 
of the search region. For this record, a dimension of 256 should be adequate for the 
current cases. But this dimension should be adjusted if other types of ducts are studied. 
A record of four of the most recent exit locations out of the lower edge should also be 


kept and checked to avoid duplicate tracking of the constant phase lines. 
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APPENDIX A: MODE LOCATIONS IN THE COMPLEX q,,/t,; PLANE 
This appendix contains figures of mode locations of evaporation ducts of 2, 4, 6, 


8, 10, 20, 30, and 40 m heights plotted in the complex q, ;/t, plane for easy comparison. 
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mode locations (2 m duct). 
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Figure A.2 qj, ;/t; mode locations (4 m duct). 
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Figure A.3 qj, ,/t, mode locations (6 m duct). 
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Figure A.4 444/44 mode locations (8 m duct) 
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Figure A.5 qj ;/t; mode locations (10 m duct). 
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(12 m duct) 
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Figure A.6 q,4/t, mode loca 
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Figure A.7 qiq/ty mode locations (14 m duct). 
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Figure A.8 qj, ,/t; mode locations (16 m duct). 
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Figure A.9 q,,/t, mode locations (18 m duct). 
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Figure A.10 qi} mode locations (20 m duct). 
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Figure A.11 qi} mode locations (22 m duct) 
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Figure A.12 q ;/t) mode locations (24 m duct). 
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Figure A.13 qj 4/t; mode locations (26 m duct). 
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Figure A.14 q,,/t; mode locations (28 m duct). 
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Figure A.15 qj ,/t; mode locations (30 m duct). 
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16 qj, ,/t; mode locations (32 m duct). 


Figure A. 
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Figure A.17 qj ty mode locations (34 m duct). 
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18 q,,/t, mode locations (36 m duct) 


Figure A. 
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Figure A.19 q), ,/t, mode locations (38 m duct). 
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Figure A.20 qj ;/t; mode locations (40 m duct). 
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APPENDIX B: SUBROUTINE FNDMOD AND FZEROX 
This appendix contains the listing of the subroutines FNDMOD and FZEROX, 
modified from the NPS version of the M-Layer FORTRAN code to track the constant 


phase-lines and to locate the modes. 
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Line# 


OPeNTWAM AWN = 


Source Line Microsoft FORTRAN Optimizing Compiler Version 5.00 


subroutine fndmod(aloss,dmmin,t1 ,dmesh,filem nrmode,geigen) 
¢ purpose. 
c This subroutine sets up the areas in the complex q11-plane 
c to search for the zeroes of the modal function. 
fe 
c inputs... 
c mxlayr - maximum number of layers allowed in refractivity 
C profile 
c nzlayr - actual number of layers in refractivity profile 
Cc aloss - maximum attenuation rate (in db/km) of modes 
c to be found 
Cc dmmin - minimum of zim(j)-zim(1) 
c tl -koal23 
c dmesh - initial mesh size divisor 
c 
c outputs... 
c qeigen - complex array containing the zeros of the modal 
c function 
c nrmode - actual number of modes found 
C 
c subroutines called... 
c fzerox 
Cc findfx 
¢ common block areas... 
Cc coml 
CREEKKEKE EE 
c 

implicit double precision (a-h,o-z) 

complex* 16 ctemp,t1,qeigen,zeros,fx1,fx2 

parameter(c20log=8.68588963806504d0,one=(1.d0,0.d0), 

$ step0=1024.d0,stepOh=step0/2.d0) 

character*3 _filem,kline 
character*40 file 

c 
Cree texters 
c¢  qeigen  - complex array containing all the zeros of the 
c modal function found 
C zeros  - complex array containing the zeros of the modal 
c function found in the current rectangular region 
c of the complex q11-plane 
c 
CREERREKE KK 
c use include file for parameters of 
c mxlayr max # layers 
Cc mxmode max # modes 
C 
Sinclude: ’mlaparm.inc’ 


***** Begin listing of: mlaparm.inc 
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“TAN WN — 


aQaaaana 


include file to define the 
maximum # of layers (mxlayr) 
maximum # of modes (mxmode) 


parameter (mxlayr=35 ) 
parameter (mxmode=127 ) 


***** End listing of: mlaparm.inc 


49 
50 
51 
3/4 


dimension geigen(mxmode),zeros(2*mxmode+ 1) 
Crete 
common /com1/waveno 
c 
ace oe ene eee ee TERE EKER EES EERE EEEE EEL ESS 
c rl - value of mode search on the real part of q!1 at the 
c left edge in tmesh units. 
c¢ 12 - value of mode search on the real part of q11 at the 
c right edge in tmesh units. 
c bot - value of the imaginary part of ql1 at the bottom 
c edge (this is set to Q). 
c 0 - value of the imaginary part of ql1 at the top edge 
Cc in tmesh units. 
c¢ step - size of search areas. 
Cee EEEKR ERE EERE AS SEE EREEEESLE LEC ERE REESE SETAE ELAVECEREAEEKA REESE 
c 
c 
CC cew ee ecew ene ce cence newer nwo ew enn nn secon enon ceeeesenesonsocosoeseces 
Cc start of executable statements 
C comme mcew ecw ewe ns cr ewnen: coc cwoew coe cnoeneccosoesecoeecesocescccsocce 
Cee TEES Eee EEE Cee eE eee Ce LEELELERELELLEELLERE EE EEE EA EREEREREERAES 
¢ set up search areas for finding modes and solve for modes, 
c calculate approximate value for re(q11) where modes start 
c occurring 
fe ee ee Ree ESAT ERT AKER EE RAE EASE A REREE ELLE EERE KERAKEE EEO ERS & 
C 
nrmode=0 
C 
recons=dreal(t1) 


rstart=-2.0d-6*dmmin 
qtest=rstart*recons 


ttop1=(2.0d-3/(waveno*c20log))*recons 
ttop=aloss*ttop1 
tmesh=ttop 1/dmesh 
if(tmesh .gt. 0.1d0) tmesh=1.0d-1 
if(tmesh .gt. 1.0d-3) then 
tol=1.0d-4 
else 
tol=tmesh*0.1d0 
end if 
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9lc¢ 


s 92 Ch**** 
93 write(*, 1002) 
94 write(16,1002) 
95 ¢ 
96 write(*,1003) tmesh,tol 
97 write(16,1003) tmesh,tol 
98 c¢ 
99 t0=dnint(ttop/tmesh)+1.d0 
100 bot=0.d0 
101 CREEK 
102 ¢c 
103 rl=dnint(qtest/tmesh)-stepOh 
104 rl=rl 
105 tr=rl+step0 
106 call findfx(rl ,t0,fx1,x1,y1,tmesh) 
107 il=int(rl) 
108 in=int(step0)+int(rl) 
109 nline=0 
110 50 continue 
111 do 100 nn=i1,in 
bi2 r2=r1+1.d0 
113 call findfx(r2,t0,fx2,x2,y2,tmesh) 
114 if (yl*y2 .lt. 0.d0) then 
1s) r=r1*tmesh 
116 write (*,5555) r 
117 c write (16,5555) r 
118 5555 _—s format (/5x,’r= ’,d28.16) 
119 go to 400 
120 endif 
121 rear 
122 fxl=fx2 
123 Xx l=xZ 
124 yl=y2 
125 100 continue 
126 retum 
127 ¢ 
128 400 continue 
129 nrold=nrmode 
130 nline=nline+1 
131 nc 100=int(nline/100) 
132 nc 10=int((nline-nc 100* 100)/10) 
133 nc l=nline-nc 100* 100-nc 10*10 
134 kline=char(48+nc 100)//char(48+nc 10)//char(48+nc 1) 
135 rfile=’e:\ans5\r’ //filem//kline//’ dat’ 
136 Crtreer 
137 open (unit=25,file=rfile,status=’unknown’) 
138 
139 call fzerox(tO,rl r11,fx1,x1,y1,fx2,x2,y2,zeros,nrold,;nmew, 
140 $ tmesh) 
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141 
142 
143 
144 
145 
146 
147 
148 
149 
150 
151 
152 
153 
154 
155 
156 
137 
158 
159 
160 
161 
162 
163 
164 
165 
166 
167 
168 
169 
170 
171 
bez 
173 
174 
i) 
176 
177 
178 
179 
180 
181 
182 
183 
184 
185 
186 
187 
188 
189 
190 


close (25) 
C*F*EEE 
Cc 

nrmode=nmew 

new0=nrmode-nrold 
Cc 
ee ee eR EEA EEE eK EREREREEE RELEASE LER SLES SE EERE LE EERE 
¢ mode search completed 
c¢ order zeros found by order of increasing real part. 
CEEEREREAERLSLELELE LER ERERERELE LER EEERERAREREREELER EEA AERERERAEEAEES 
Cc 

if(nrmode .gt. 1) then 

jkend=nrmode-1 


do 420 jk=1 jkend 
nend=nrmode-jk 


do 410 ja=1,nend 
nrj=nrmode-ja 
nrjl=nrj+1 
if(dimag(zeros(nrj1)).It. 
$ dimag(zeros(nrj))) then 
ctemp=zeros(nrj1) 
zeros(nrj1)=zeros(nrj) 
zeros(nrj)=ctemp 
end if 
410 continue 
420 _ continue 
end if 
Cc 
CRRRREREKK EER ERR EKEK KKK KKK EK KEEKK KEK KK KEE EK ERK KK KK KEKE EKKEKKKEKEKKKEKE KK KX 
c_ the possibility exists that duplicate (within the tolerance ‘tol’) 
c zeros of the modal equation will be found. eliminate these 
c duplicate zeros. 
CERRRREKAAE ERK EKER EK KEKE KK KK KEK KAK KKK KA EERE KKK KK KEKE KKK KKK RE K KK KKK K KKK 
Cc 


jkflag=0 
jkend=nrmode-1 


do 240 jk=1,jkend 
jkl=jk+1 
ctemp=zeros(jk1) 
chksq=cdabs(zeros(jk-jkflag)-ctemp) 
if(chksq .lt. tol) then 
jkflag=jkflag+1 
go to 240 
end if 
zeros(jk 1-jkflag)=ctemp 
240 continue 
nrmode=nrmode-jkflag 
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191 
192 
193 
194 
195 
196 
197 
198 
19 


201 
202 
203 
204 
205 
206 
207 
208 
209 
210 
211 
212 
213 
214 
215 
216 
217 
218 
219 
220 
221 
@22 
223 
224 
225 
226 
Zo) 
228 
229 
230 
Zo 
232 
259 
234 
Zo) 
236 
25] 
238 
239 
240 


Cc 
CREERERESEERS ERSTE EEEEEE EEE EESEL EEE NES EEO 


c Store the zeros as the eigenvalues. 
CUEEREAREL ADEA LS CHESS SES KEES ES Re 


c 
nrmode=min0(nrmode,mxmode) 
do 430 jk=1,nrmode 
qeigen(jk)=zeros(jk) 
430 continue 


Cc 

il=int(rl1)+1 

if (il .It. in) then 

rl=rl1+1.d0 

call findfx(r1,t0,fx1,x1,y1,tmesh) 

go to S50 

endif 

return 
Cc 
Cc eceecccees cor eree ee ece ce ee ee ee cesses ee sesccos seesesssesesseosecseesecs 
Cc format statements 
Cc eecececwccccco cre swocecowecerececssesacosenscsesscosesessssosocesssscco 
1000 format(/Sx,’searching for zeroes in this areas are’, 


$ ° defined by:’/Sx,’istart= ’,i10/5x,’itop= ’,i10, 
$ 5x,” ibot= ’,111/5x,’ileft= °,110,5x,’iright= ’,110) 


1001 format(5x,i4,’ new zeroes found in this area.’/) 


1002 format(//Sx,’'******* start search for modal eigenvalues’, 


$ 2. eee ee) 


1003 format(/Sx,’tmesh= ’,d15.5,5x,’tol= ’,d15.5/) 


end 
Cc 
Cc 
CERRRERREKEREKEKEEREER EKER EEEREEKEKEEEEEKEEREKEREREREEREKEERE EERE AEEREES 


CRRERREREERER EERE EERE A RER HERE RERKERREK EEK ERERRERER EEE SEE EEE EEEKKKEEE SE 


C 
subroutine fzerox(t0,rl,rl 1,fx1,x 1,y1,fx2,x2,y2,zeros,nrold, 
$ nmew,tmsh0) 
C 
C#*EE* 
c fzerox is a routine for finding the zeroes of a complex function, f, 
c¢ which lie within a specified rectangular region of the 
complex ql! plane, assuming that the function has only 
simple zeroes over this rectangle. 


parameters specifying the search rectangle: 


Cc 
C 
c 
Cc 
¢  tmesh - set equal to about half the average spacing between 
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241 
242 
243 
244 
245 
246 
247 
248 
249 
250 
251 
252 
253 
254 
255 
256 


zeroes within the rectangle. A smaller value may be used 
as a safety measure, but too small a value will result 
in excessively long run time. 
zeros - Output list of (complex) values of ql1 at which 
zeroes are found. 
nmew,nrold- the number of zeroes found 


subroutines called-- 
findfx 
roots 


Qa aQaaaananaaannn 


Cree 
implicit double precision (a-h,o-z) 
complex*16 £10,f01,f11,fx 1,fx2,fx00,fx 10,fx01,fx11, 
+ one,sol,zeros 
parameter(one=(1.d0,0.d0)) 
Sinclude: *mlaparm.inc’ 


***** Begin listing of: mlaparm.inc 


“ANNA RP WY — 


Cc 

Cc include file to define the 

c maximum # of layers (mxlayr) 

Cc maximum # of modes (mxmode) 
Cc 


parameter (mxlayr=35 ) 
parameter (mxmode=127 ) 


**%** End listing of: mlaparm.inc 


paud| 
258 
259 
260 
261 
262 
263 
264 
265 
266 
267 
268 
269 
270 
Fai | 
272 
275 
274 
295 
276 
Ped | 
278 
279 
280 
281 


dimension sol(3),theta(2),zeros(2*mxmode+ 1),rline(2,1024) 
c 
common/tmccom/tmesh 
Ce ERE 
Cc 
npoint=0 
nmew=nrold 
tmesh=tmsh0 
bot=0.d0 
t=t0-1.d0 
rll=rl 
rrll 
C*F*E* 
Cc 
fxOl=fx1 
fxOlr=x1 
fxO0li=y1 
ixl=txZ 
fx1lr=x2 
fxlli=y2 
go to 82 
Cc 
C#RE*E* 
c enter mesh square from left side or exit rectangle at right edge. 
c 
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331 


20 ~=sr=r+1.d0 

21 = fx01l=fx11 
fxOlr=fx11r 
fx01i=fx1 li 
fx00=fx 10 
fxO0Or=fx 10r 
fxO0i=fx10i 

22 ~=continue 


call findfx(r+1.d0,t+1.d0,fx11,fx 1 1r,fx11i,tmesh) 
call findfx(r+1.d0,t.fx 10,fx 10r,fx10i,tmesh) 
Cette 
c Determine the edge of exit of im(f)=0 from current mesh. 


edgeit=fx01i*fx11i 
edgeib=fx00i* fx 10i 
if (edgeib .gt. 0.d0) then 
e Im(f)=0 goes through the 01 to 10 line. 
if (edgeit .gt. 0.dO) then 
c Im(f)=0 goes through the 10 to 11 edge (edge 1). 
lout=1 
else 
Cc Im(f)=0 goes through the 01 to 11 edge (edge 2) 
lout=2 
end if 
else 
c Im(f)=0 goes through the 00 to 10 edge (edge 4) 
lout=4 
if (edgeit .It. 0.d0) then 
c Im(f)=0 also runs through 01 to 11 and 10 to 11 edges. 


C Store crossing location and in/out information. 
knot34=knot34+1 
C loc34r(knot34)=lIr 
Cc loc34i(knot34)=li 
end if 
end if 
CEREREES 
go to 85 
Ctreet 
c enter mesh square from bottom side or exit rectangle at top edge. 


40 t=t+1.d0 
41 fx00=fx01 
fxO0Or=fx01r 
fx00i=fx01i 
fx10=fx11 
fxl0r=fx11r 
fx10i=fx11i 
42 continue 
call findfx(r,t+1.d0,fx01.fx01r,fx01i,tmesh) 
call findfx(r+1.d0,t+1.d0,fx11,fx 1 1r,fx11i,tmesh) 


Ctr rene 


c¢ Determine the edge of exit of im(f)=0 from current mesh. 
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edgeil=fx00i* fx01i 
edgeir=fx 10i*fx11i 
if (edgeir .gt. 0.dO) then 
Im(f)=0 goes through the 00 to 11 line. 
if (edgeil .gt. 0.dO) then 
Im(f)=0 goes through the 01 to 11 edge (edge 2) 
lout=2 
else 
Im(f)=0 goes through the 00 to 01 edge (edge 3). 
lout=3 
end if 
else 
Im(f)=0 goes through the 10 to 11 edge (edge 1) 
lout=1 
if (edgeil .lt. 0.dO) then 
Im(f)=0 also runs through 00 to 01 and 01 to 11 edges. 
Store crossing location and in/out information. 
knot4 l=knot41+1 
loc4 Ir(knot41)=Ir 
loc4 li(knot4 1)=li 
end if 
end if 


CRReeeee 


go to 85 


errnre 


¢ enter mesh square from right side or exit rectangle at left edge. 


60 = r=r-1.d0 


61 
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fx11l=fx01 
fxllr=fxOlr 
fx1 li=fx01i 
fx 10=fx00 
fx 1Or=fx00r 
fx10i=fx00i 
continue 
call findfx(r,t+1.d0,fx01,fx01r,fx01i,tmesh) 
call findfx(r,t,fx00,fx00r,fx00i,tmesh) 


Cut eeeee 


c 


Determine the edge of exit of im(f)=0 from current mesh. 
edgeit=fx01i*fx1li 
edgeib=fx00i* fx 10i 
if (edgeit .gt. 0.dO) then 
Im(f)=0 goes through the 01 to 10 line. 
if (edgeib .gt. 0.d0) then 
Im(f)=0 goes through the 00 to 01 edge (edge 3). 
lout=3 
else 
Im(f)=0 goes through the 00 to 10 edge (edge 4) 
lout=4 
end if 


Sf 


else 
c Im(f)=0 goes through the 01 to 11 edge (edge 2) 
lout=2 
if (edgeib .It. 0.dO) then 
é Im(f)=0 also runs through 00 to 10 and 00 to O1 edges. 
(© Store crossing location and in/out information. 
knot12=knot12+1 
Cc loc12r(knot12)=Ir 
C loc 12i(knot12)=hi 
end if 
end if 
ChRRRKAE 
go to 85 
CXFER* 
c enter mesh square top side or exit rectangle at bottom edge. 


c 

80 t=t-1.d0 

81  fx0l=fx00 
fx0 Ir=fx00r 
fx0 1li=fx00i 
fxll=fx10 
fxl lr=fx10r 
fx1 li=fx10i 

82 continue 


call findfx(r,t,£x00,fx00r,fx001,tmesh) 
call findfx(r+1.d0,t,fx 10,fx 10r,fx 10i,tmesh) 
c 
crtteree 
c Determine the edge of exit of im(f)=0 from current mesh. 


C 
83 continue 
edgeil=fx00i*fx01i 
edgeir=fx 10i*fx11i 
if (edgeil .gt. 0.dO) then 
Cc Im(f)=0 goes through the 00 to 11 line. 
if (edgeir .gt. 0.d0) then 
Cc Im(f)=0 goes through the 00 to 10 edge (edge 4) 
lout=4 
else 
Cc Im(f)=0 goes through the 10 to 11 edge (edge 1). 
lout=1 
end if 
else 
‘s Im(f)=0 goes through the 00 to 01 edge (edge 3) 
lout=3 
if (edgeir .It. 0.d0) then 
Cc Im(f)=0 also runs through 00 to 10 and 10 to 11 edges. 
Store crossing location and in/out information. 
knot23=knot23+1 
fe loc23r(knot23)=Ir 


Q 
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461 
462 
463 
464 
465 
466 
467 
468 
469 
470 
471 
472 
473 
474 
475 
476 
477 
478 
479 
480 
481 


Cc loc23i(knot23)=li 
end if 
end if 


CERREREAARAEKREE AH EREREA EERE EEE LEER ELE ERE LEE EAE HAR AAAS EEE ER EEE ERLE EEE EE 


c Test for there being at least one re(f)=0 line entering and 


c leaving the mesh square. 
CERRRAEAEARRAAKAA HRA RRAAAER RELEASE BEER LHR EER AAA AKER A HAE HARRAH EERE EE EEE 
c 
85 continue 
npoint=npoint+ 1 
rir=r*tmesh 
ttt=t*tmesh 
rline(1 npoint)=rr 
rline(2,npoint)=ttt 
if (npoint .eq. 1024) lout=5 
Cc 
if(t .It. bot) .or. (t .gt. t0)) go to 100 
if ((fxOOr*fx10r .gt. 0.d0) .and. (fxO1r*fx1 Ir .gt. 0.d0) 
+ .and. (fx00r*fxOlr .gt. 0.d0)) go to (20,40,60,80, 100) 
+ lout 
CEE EERER KEEEEE CLES EEE ECELELELELE MELEE KEKE RAR AEA AAAL LAA LE RELLELE RAK KE HAS 


c Computate the values of the modal function at the corners of 
¢ amesh square to determine its Taylor series to the 3rd order 


c for estimating its root locations. 
ee ee ee OO ee CE EEEERE ELE EEE ESE ELE RE EAEAL KE LELERKERAEEE ES & 
Cc 
c  f00=one 
f10=cdexp(fx 10-fx00)-one 
f01=cdexp(fx01-fx00)-one 
f1 1=cdexp(fx11-fx00)-one 
c 
Cutt ttesere* estimate locations of zeroes by radicals***************** 


C 
call roots(f10,f01,f11,sol,nrsol) 


do 90 n=1,nrsol 
ureal = dreal(sol(n)) 
uimag = dimag(sol(n)) 
if (ureal .It. 0.dO .or. ureal .gt. 1.0d0) go to 90 
if (uimag It. 0.dO .or. uimag .gt. 1.0d0) go to 90 
92 theta(1)=(r+ureal)*tmesh 
theta(2)=(t+uimag)*tmesh 
nmew = nmew+1 
zeros(nrnew)=dcmpIx(theta(1),theta(2)) 
90 continue 
Cc 
CRE RHHRH H H A H HHH AR HR HR HH HH AH AH 


¢ continue following the phase line 
C8 HHH AH HA HH HH RH HR HH HH HH HH A 
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482 c 

483 95 continue 

484 go to (20,40,60,80,100) lout 

235 c** 

486 100 continue 

487 Chtteee 

488 write (25,8888) ((rline(i,j), i=1, 2), j=1, npoint) 
489 8888 format (2e15.5) 


490 c 
491 retum 
492 end 


60 


APPENDIX C: CONSTANT PHASE LINE TRACKING FLOW CHARTS 


nrmode = 0 

recons = dreal(t1) 

qtest=-2x10*(-6) x dmmin x recons 

ttop1=((2x10*(-3) / k x 20log(e)) x 
recons 

top = aloss x ttop1 

tmesh = ttop1 / dmesh 


tmesh = 0.1 
tS ore 
tol = tmesh x 0.1 tol = 0.0001 


(tmesh, tol) 
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t0=dnint(ttop/tmesh)+1.dO 
bot = 0.d0 


r1=dnint(qtest/tmesh)-stepOh 
rl=r1 
rr=r1+1024.d0 


CALL FINDFX 


11 = int(r1) ,in = 1024+int(r1), nline=0 


DO 100 
nn=!1,in 
r2=r1+1.d0 


N M4 
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nrold = nrmode 

nline = nline + 1 

nc100 = int(nline/100) 

neiO = int((nline-nc100*100)/10) 

nei =niine - nc100*100 - nc10*10 
kline=char(48+nc100)//char(48+nc10 

)/[char(48+nc1) 

rfile='e:\mi.che\ans\r'//filem//kline//.dat' 


Open(unit=25, file=rfile, status= 
‘unknown’') 


CALL FZEROX 


CLOSE (25) 


nrmode =nrnew 
new0 = nrmode - nrold 
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WRITE 
( nrnew ) 


DO 420 
jk =1, jkend 
nend = nrmode - jk 


DO 410 

ja= 1, nend 
nj = nrmode -ja 
ni =nry+1 


dreal(zero(nrj1) 
<dreal(zeros(nr)) 


ctemp = zeros (ni1) 
zeros (nrj1) = zeros (nr) 
zeros (nr) = ctemp 
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kflag=jkflag+1 


CONTINUE 


DO 240 

Jk = 1, jkond 
jkiejk+1 
ctamp = zeros (jk1) 


chksq = 
cdabs(zeros 
(jk-jkflag)- 
ctemp) 


zaros()k1-jkflag) = ctamp 


nrmode=nrmode - 
jkflag 
nrmode=minO(nrmode, 
mxmode) 


DO 430 
jk = 1, nrmode 


qeigen(jk) = zeros(jk 
it = int(r11) + 1 
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<i ) 
Y 
ri =r11 + 1.d0 
CALL FINDFX 
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APPENDIX D: CONSTANT PHASE LINES IN THE COMPLEX q, ,/t,; PLANE 
This appendix contains plots of the constant phase lines Im{D(q)} =0 and 
Re{D(q)} = 0 initiating from the top edge of the search region in the complex q,4/ty 


plane for the evaporation ducts of 2, 4, 6, 8, 10, 20, 30, and 40 m heights. 
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Figure D.1 Constant phase lines in the q, ,/t, plane (2 m duct). 
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the q, ,/t, plane (4 m duct) 
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Figure D.2 Constant phase 
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D.3 Constant phase lines in the q, ;/t, plane (6 m duct). 


Figure 
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Figure D.4 Constant phase lines in the qj /ty plane (8 m duct). 
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the q14/t; plane (10 m duct) 


lines in 


Figure D.5 Constant phase 





Figure D.6 Constant phase lines in the q,/t, plane (20 m duct). 
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Figure D.7 Constant phase lines in the 44 4/t; plane (30 m duct). 
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Figure D.8 Constant phase lines in the q11/t; plane (40 m duct). 
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APPENDIX E: SEPARATION OF Im{D(q)}=0 LINES 
This appendix contains figures of the separation of Im{D(q)}=0 lines along the top 
edge of the search region. The 2, 4, 6, 8, 10, 20, 30 and 40 meter evaporation ducts are 


included. 
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Figure E.1 Separation of Im{D(q)}=0 lines along the top edge of the search region 
(2 m duct). 
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Figure E.2 Separation of Im{D(q)}=0 lines along the top edge of the search region 
(4 m duct). 
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Figure E.3 Separation of Im{D(q)}=0 lines along the top edge of the search region 
(6 m duct). 
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Figure E.4 Separation of Im{D(q)}=0 lines along the top edge of the search region 
(8 m duct). 


79 


(1/1 [bay Jo sduarajsIG 


Figure E.5 Separation of Im{D(@q)}=0 lines along the top edge of the search region 
(10 m duct). 
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Figure E.6 Separation of Im{D(q)}=0 lines along the top edge of the search region 
(20 m duct). 
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Figure E.7 Separation of Im{D(q)}=0 lines along the top edge of the search region 
(30 m duct). 
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Figure E.8 Separation of Im{D(q)}=0 lines along the top edge of the search region 
(40 m duct). 
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APPENDIX F: DIFFERENCE BETWEEN CONSECUTIVE Re(q,,/t;) VALUES 
Separation in real part of neighboring q-eigenvalues scaled by t, is plotted in this 


appendix against the eigenvalue locations along the real q, ;/t; axis. 
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Figure F.1 Difference between consecutive Re(q,/t)) values (2 m duct). 
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Figure F.2 Difference between consecutive Re(q,,/t)) values (4 m duct). 
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Figure F.3 Difference between consecutive Re(q,,/t)) values (6 m duct). 
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Re{qm/tl } (6 m duct) 
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Figure F.4 Difference between consecutive Re(q,/t,) values (8 m duct). 


88 


ss 
=! 
a) 
5 
sc 
= 
© 
=— 
— 
2 
oy 
o 
a 


a) 


a @e@ © t+ A © 
mt mm =o 


= — = SS 
Ss So 5 


anjeA ([i/wb)ay sannoasuod usemjog souaIIIJIG 





Figure F.5 Difference between consecutive Re(q,/t) values (10 m duct). 
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Figure F.6 Difference between consecutive Re(q,,/t)) values (20 m duct). 
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Figure F.7 Difference between consecutive Re(q,/t)) values (30 m duct). 
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Figure F.8 Difference between consecutive Re(q,,/t)) values (40 m duct). 
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